#!/bin/bash
set -e
# no reference
# only support window=20,000 window=200,000

while getopts  ":w:c:p:b:f:" opts
do
	case  $opts  in
	w)
		window=$OPTARG
		;;
	c) chrom=$OPTARG;;
	b)
        asemmbly=$OPTARG
        ;;
	p)
        out_prefix=$OPTARG
        ;;
	f)
        filt=true
		f_arg=$OPTARG
        ;;
	\?)
		echo `basename $0` "[-w window_size] [-b asemmbly] [-p out_prefix] [-f filter_window_size] <sample.bam>"
		exit 1
		;;
	esac
done
shift $(($OPTIND - 1))

if [ -z "$1" ]; then
	echo `basename $0` '[-w window_size] [-b asemmbly] [-p out_prefix] [-f filter_window_size] <sample.bam>'
	exit 1
fi


#-----------------------------------------------
. /mnt/ilustre/app/medical/tools/.var #---------
#-----------------------------------------------


if [ -z "$window" ]; then
	window=200000
fi

if [ -z "$f_arg" ]; then
	f_arg=1
fi


if
samtools view -H $1 |grep '^@RG'
then
sample_name=`samtools view -H $1 |grep '^@RG' |perl -n -e 'if(/SM:(.+?)[\t|\n]/){print $1}'`
fi
echo $sample_name

if test -z $chrom; then
chrom=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,X,Y
fi
echo $chrom

base_ref_genome=`basename $ref_genome`

if [ -e "${ref_genome}.gc.${window}.wig" ]; then
	echo
	echo ${ref_genome}.gc.${window}.wig exists.
else
	echo
	echo
	echo gcCounter
	echo window $window
	echo $ref_genome
	gcCounter -c $chrom -w $window $ref_genome > $ref_genome.gc.$window.wig
fi

if [ -e "${ref_genome}.map100bp.${window}.wig" ]; then
	echo ${ref_genome}.map100bp.${window}.wig exist.
else
	echo
	echo
	echo mapCounter
	mapCounter -c $chrom -w $window ${ref_genome}.map100bp.bw \
	> ${ref_genome}.map100bp.${window}.wig
fi

bai1=$1.bai
# bai2=$2.bai
base1=`basename $1 .bam`
# base2=`basename $2 .bam`
# control_dir=`dirname $2`
if [ ! -e "$bai1" ]; then
	
	echo samtools index $1
	samtools index $1
fi


if test ! -e $out_prefix.$window.sample.wig
then
echo
echo
echo readCounter sample
echo window $window
echo $1
readCounter -c $chrom -w $window $1 > $out_prefix.$window.sample.wig
else 
echo
echo
echo $out_prefix.$window.sample.wig exists.
fi


Rscript ${tools_path}/script/hmmcopy0901.r --args \
$out_prefix.$window.sample.wig \
$sample_name \
$ref_genome.gc.$window.wig \
$ref_genome.map100bp.$window.wig \
$chrom

mv 1.cnv.hmmcopy.pdf $out_prefix.cnv.hmmcopy.$window.pdf

if test "$filt" = true; then

#-----------------------
cp ${tools_path}/matlab/rc_filt1.m .

cat $out_prefix.$window.sample.wig |awk '{if(/^f/){print "-1"}else{print}}' > 1.rc.txt

echo;echo;echo matlab should run foreground'!'
matlab -nojvm -r "rc_filt1; quit;"

cat 1.rc_filt1.txt |awk 'BEGIN{chr="'$chrom'"; split(chr, chr1, ","); i=0}{if(/^-/){i++; print "fixedStep chrom=" chr1[i] " start=1 step=200000 span=200000"}else{print}}' > $out_prefix.$window.sample.filt1.wig

# paste $out_prefix.$window.sample.wig 1.rc_f.txt > 1.merge.wig.txt

# cat 1.merge.wig.txt |awk -F "\t" '{if(/^f/){print $1}else{print $2}}' > $out_prefix.$window.sample.filt.wig

#-----------------------

window1=$[window * 10]

sample_name1=${sample_name}_f$f_arg

Rscript ${tools_path}/script/hmmcopy0901.r --args \
$out_prefix.$window.sample.filt1.wig \
$sample_name1 \
$ref_genome.gc.$window1.wig \
$ref_genome.map100bp.${window1}.wig \
$chrom

mv 1.cnv.hmmcopy.pdf $out_prefix.cnv.hmmcopy.$window.$f_arg.pdf

fi
